Rasch analysis of outcomes

PDSS-SR

Authors
Affiliations
Nils Hentati Isacsson

Centre for Psychiatry Research, Department of Clinical Neuroscience, Karolinska Institutet, & Stockholm Health Care Services, Region Stockholm, Sweden

Magnus Johansson
Published

2024-12-19

1 Setup

Code
library(easyRasch) # devtools::install_github("pgmj/easyRasch")
#library(RISEkbmRasch)
library(grateful)
library(ggrepel)
library(car)
library(kableExtra)
library(tidyverse)
library(eRm)
library(iarm)
library(mirt)
library(psych)
library(psychotree)
library(matrixStats)
library(reshape)
library(knitr)
library(patchwork)
library(formattable) 
library(glue)


### optional libraries
#library(TAM)
#library(skimr)
#library(janitor)
# 1727699
# nohup quarto render PDSSSR.qmd --to html &> render_pdsssr.out & 

### some commands exist in multiple packages, here we define preferred ones that are frequently used
select <- dplyr::select
count <- dplyr::count
recode <- car::recode
rename <- dplyr::rename

path_prim <- '/srv/projects/nils/study4rasch'
path2 <- '/Volumes/projects/k8_CPF_Kaldo/Data/Lärande Maskiner/Nils/Study4data'
running_machine = Sys.info()[['sysname']]
path = ifelse(running_machine!='Linux',path2,path_prim)
set_mywd_path = ifelse(running_machine!='Linux','./study4rasch/rasch_measures','./rasch_measures')
n_cores =ifelse(running_machine!='Linux',1,8)

tryCatch({setwd(set_mywd_path)}, error = function(set_mywd_path){print('already in right directory')})
[1] "already in right directory"
Code
source('settings.R') # containing item labels 
source('mod_easyRasch_func.R') # modified functions + added 

2 Importing data

Code
### import data
df_raw <- read.csv(file=file.path(path,'data','raw_data','outcome.csv')) #

### Clean data 
source("preprocess.R")
df <- preprocess_data(df_raw)

# items to use 
items_to_use <- 'PDSS.SR' # options are #MADRS.S, LSAS.SR, PDSS.SR
itemlabels <- itemlabels %>% 
  as_tibble() %>%
  filter(str_detect(itemnr, items_to_use))


# DIF+aux variables to use 
# DIF not used 'marital_status','children', 'working'
dif_variables <- c('Patient','Treatment','sex','age','TreatmentAccessStart','education')

### Make a backup of the dataframe, in case you need to revert changes at some point
df.all <- df

print(itemlabels)
# A tibble: 7 × 2
  itemnr    item                
  <chr>     <chr>               
1 PDSS.SR_1 Number              
2 PDSS.SR_2 Distress            
3 PDSS.SR_3 Worry/Fear          
4 PDSS.SR_4 Avoidance places    
5 PDSS.SR_5 Avoidance activity  
6 PDSS.SR_6 Impairment work/home
7 PDSS.SR_7 Impairment social   

3 Checking missing

Code
itemStart <- items_to_use
out = RImissing_mod(df,itemStart,facet="Time")
print(out)

As we see here PDSS.SR is not used for all conditions over the course of their treatments. However during screening they are. We start by showing all timepoints, to draw targeting plots.

3.1 Remove missing

Code
##### Before filtering out participants, you should check the missing data structure using RImissing() and RImissingP()

# If you want to include participants with missing data, input the minimum number of items responses that a participant should have to be included in the analysis:
min.responses <- 4 # In our case if they have missing on one item all items for that time are missing. 

df_save <- df
# Select the variables we will work with, and filter out respondents with missing data
df <- df %>% 
  select(all_of(c(dif_variables,"Time",itemlabels$itemnr))) %>%  # v
  filter(rowSums(is.na(select(.,all_of(itemlabels$itemnr)))) < min.responses) 

4 Create DIF variables

Code
#---- Create DIF variables----
  
# DIF variables into vectors, recoded as factors since DIF functions need this
# these could also be stored in its own dataframe (not a tibble) instead of as vectors
# Named vector for the new types

type_transform <- c(Treatment = "factor", sex = "factor", age = "numeric",
                    TreatmentAccessStart ="integer",education="factor",Time='factor')

# Transform columns based on the named vector
dif <- df %>%
  mutate(across(names(type_transform), ~ switch(type_transform[which(names(type_transform) == cur_column())], 
                                               "integer" = as.integer(.), 
                                               "character" = as.character(.),
                                               "factor" = as.factor(.), 
                                               "numeric" = as.numeric(.)))) %>%
                                               as.data.frame(.) %>%
                                               select(!all_of(itemlabels$itemnr))


# then remove them from dataframe, since we need a dataframe with only item data for the Rasch analyses
df <- df %>% select(all_of(c("Time",itemlabels$itemnr))) # add time here ? 
dfnotime <- df %>% select(!Time)
source("RISE_theme.R")

5 Demographics descriptives

Code
dif_spec <- dif %>% filter(Time=='SCREEN') %>% select(!all_of(c("Time","Patient")))
summary(dif_spec)
# RIdemographics(dif_spec,'Treatment') <- function crashes TODO make own 
 Treatment  sex           age        TreatmentAccessStart         education   
 MDD:2853   F:3834   Min.   :18.00   Min.   :2008         Primary      : 484  
 PD :1592   M:2366   1st Qu.:27.00   1st Qu.:2012         Secondary    :3024  
 SAD:1755            Median :33.00   Median :2014         Postsecondary:2692  
                     Mean   :35.24   Mean   :2014                             
                     3rd Qu.:42.00   3rd Qu.:2017                             
                     Max.   :84.00   Max.   :2019                             

6 Overall number of responses

Code
# Collapsed 
RIallresp(dfnotime)
# Seperate 
RIallresp_over_times <- df %>% # 
  split(.$Time) %>% # split the data 
  map(~ RIallresp(.x %>% dplyr::select(!Time))) #+ labs(title = .x$Time)) # create separate for each time

# make nice later 
RI_allresp_kable_grid = combine_kables_grid(RIallresp_over_times,cols=3)
RI_allresp_kable_grid
Response category
Number of responses
Percent
0
48333
31.8
1
54744
36.1
2
34135
22.5
3
12359
8.1
4
2182
1.4
x
SCREEN.Response category SCREEN.Number of responses SCREEN.Percent PRE.Response category PRE.Number of responses PRE.Percent WEEK01.Response category WEEK01.Number of responses WEEK01.Percent
0 16055 37.0 0 1661 14.1 0 1915 19.8
1 10996 25.3 1 3725 31.7 1 3879 40.1
2 10206 23.5 2 4114 35.0 2 2868 29.6
3 5057 11.7 3 1921 16.4 3 904 9.3
4 1086 2.5 4 325 2.8 4 115 1.2
WEEK02.Response category WEEK02.Number of responses WEEK02.Percent WEEK03.Response category WEEK03.Number of responses WEEK03.Percent WEEK04.Response category WEEK04.Number of responses WEEK04.Percent
0 2114 21.6 0 2413 24.7 0 2649 28.0
1 4268 43.6 1 4262 43.6 1 4041 42.8
2 2545 26.0 2 2355 24.1 2 2093 22.1
3 768 7.8 3 665 6.8 3 573 6.1
4 98 1.0 4 70 0.7 4 94 1.0
WEEK05.Response category WEEK05.Number of responses WEEK05.Percent WEEK06.Response category WEEK06.Number of responses WEEK06.Percent WEEK07.Response category WEEK07.Number of responses WEEK07.Percent
0 2679 29.6 0 2903 33.3 0 2841 34.1
1 3911 43.2 1 3656 41.9 1 3482 41.8
2 1887 20.8 2 1669 19.1 2 1540 18.5
3 487 5.4 3 425 4.9 3 397 4.8
4 87 1.0 4 76 0.9 4 70 0.8
WEEK08.Response category WEEK08.Number of responses WEEK08.Percent WEEK09.Response category WEEK09.Number of responses WEEK09.Percent WEEK10.Response category WEEK10.Number of responses WEEK10.Percent
0 3020 37.8 0 2962 38.8 0 3027 42.3
1 3273 40.9 1 3082 40.4 1 2842 39.7
2 1348 16.9 2 1230 16.1 2 1042 14.6
3 313 3.9 3 325 4.3 3 220 3.1
4 40 0.5 4 31 0.4 4 30 0.4
POST.Response category POST.Number of responses POST.Percent " " " "
0 4094 45.4
1 3327 36.9
2 1238 13.7
3 304 3.4
4 60 0.7

7 Descriptives of raw data

Response distribution for all items are summarized below.

7.1 Floor/ceiling effects

Code
RIrawdist(dfnotime)

Code
#df_test = df %>% filter(Time=='SCREEN') %>% select(!Time)
# Across time 

# seperated time - might need to fix this remake in ggplot?
png(filename=paste0("./plots/",itemStart,"_rawdist_over_time.png"),
width=800,height=2400)
par(mfrow=c(7,2),mar = c(4, 4, 4, 2))  # c(5 <- bottom lines, 4, 4, 2)
RIrawdist_over_times <- df %>% # 
  split(.$Time) %>% # split the data 
  imap(~ {
    RIrawdist(.x %>% dplyr::select(!Time))
    mtext(paste("Time:", .y), side = 1, line = 3, adj = 0, cex = 0.8)
    })

dev.off()
png 
  2 
Code
print('Saved plot to /plots')
[1] "Saved plot to /plots"

7.2 Guttman structure

Code
RIheatmap(dfnotime)

Code
ggplots_with_spec_func_over_time(df,func_call = RIheatmap,title_all='Guttman structure over time',blankx=TRUE)

7.3 Descriptives - item level

Code
df_test = df %>% filter(Time=='SCREEN') %>% select(!Time)
RIlistItemsMargin(df, fontsize = 12)
itemnr item
PDSS.SR_1 Number
PDSS.SR_2 Distress
PDSS.SR_3 Worry/Fear
PDSS.SR_4 Avoidance places
PDSS.SR_5 Avoidance activity
PDSS.SR_6 Impairment work/home
PDSS.SR_7 Impairment social
Code
# over all times 
RItileplot(dfnotime)

Code
# by time 
ggplots_with_spec_func_over_time(df,func_call = RItileplot,title_all='',blankx=FALSE)

Code
RIbarstack(dfnotime)

Code
# by time 
ggplots_with_spec_func_over_time(df,func_call = RIbarstack,title_all='',blankx=FALSE)

Code
RIbarplot(dfnotime)

8 Rasch analysis 1

Code
df_test = df %>% filter(Time=='SCREEN') %>% select(!Time)
RIlistItemsMargin(df, fontsize = 12)
itemnr item
PDSS.SR_1 Number
PDSS.SR_2 Distress
PDSS.SR_3 Worry/Fear
PDSS.SR_4 Avoidance places
PDSS.SR_5 Avoidance activity
PDSS.SR_6 Impairment work/home
PDSS.SR_7 Impairment social

The eRm package, which uses Conditional Maximum Likelihood (CML) estimation, will be used primarily. For this analysis, the Partial Credit Model will be used.

8.1 Timepoint decision

Due to the amount of data across timepoints targeting will first be inspected, and the best fitting timepoint will be chosen based on targeting, which in turn will inform the rest of the analyses.

8.1.1 Targeting

Code
# increase fig-height above as needed, if you have many items og. value was 5. 
#RItargeting(dfnotime)

require(patchwork)
targeting_time_plots <- ggplots_with_spec_func_over_time(
  df,func_call = RItargeting,title_all='',
  blankx=FALSE,colsplot = 3,rowsplot=5,return_list_plots=TRUE)
good_layout_targeting_plots <- wrap_plots(targeting_time_plots, ncol = 3) + plot_layout(guides= "collect")
#print(good_layout_targeting_plots)
ggsave(
  filename = paste0("./plots/",items_to_use,"_targeting_over_time.png"),
  plot = good_layout_targeting_plots,
  width = 20,  # Increase width (in inches)
  height = 45, # Increase height (in inches)
  dpi = 300    # High resolution
)
print(good_layout_targeting_plots)

Pre has very good targeting looks promising. We also note that this is only for patient in treatment for panic disorder (no other treatment-conditions are mixed in PRE). We doublecheck how much data we have for each threshold. (Recomm. 20+)

8.1.2 Number of responses per threshhold

Code
dataframen4 <- df %>% filter(Time=="PRE") %>% select(!Time)
sapply(dataframen4, function(column) table(column))
  PDSS.SR_1 PDSS.SR_2 PDSS.SR_3 PDSS.SR_4 PDSS.SR_5 PDSS.SR_6 PDSS.SR_7
0       129       157        52       161       521       368       273
1       673       389       409       504       526       621       603
2       679       807       628       660       430       426       484
3       148       296       496       293       178       227       283
4        49        29        93        60        23        36        35

Based on the above we move forward with the PRE timepoint.

8.2 Analysis continuation PRE

Code
df_pre <- df %>% filter(Time=='PRE') %>% select(!Time) #%>% sample_n(800)
Code
RIitemfit(df_pre, cutoff = "Smith98")
Item InfitMSQ Location 1 ± 2 / √n
PDSS.SR_1 0.996 -0.10 [0.951, 1.049]
PDSS.SR_2 0.875 -0.07 [0.951, 1.049]
PDSS.SR_3 0.973 -0.96 [0.951, 1.049]
PDSS.SR_4 1.045 -0.24 [0.951, 1.049]
PDSS.SR_5 1.29 0.75 [0.951, 1.049]
PDSS.SR_6 0.956 0.40 [0.951, 1.049]
PDSS.SR_7 0.915 0.22 [0.951, 1.049]
Note:
MSQ values based on conditional estimation (n = 1678 complete cases).
Code
simfit1 <- RIgetfit(df_pre, iterations = 300, cpu = n_cores)  

RIitemfit(df_pre, simfit1)
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Location
PDSS.SR_1 0.996 [0.929, 1.081] 0.981 [0.927, 1.075] no misfit no misfit -0.10
PDSS.SR_2 0.875 [0.921, 1.085] 0.874 [0.918, 1.091] 0.046 0.044 -0.07
PDSS.SR_3 0.973 [0.946, 1.06] 0.976 [0.946, 1.065] no misfit no misfit -0.96
PDSS.SR_4 1.045 [0.919, 1.084] 1.057 [0.914, 1.088] no misfit no misfit -0.24
PDSS.SR_5 1.29 [0.899, 1.093] 1.342 [0.875, 1.131] 0.197 0.211 0.75
PDSS.SR_6 0.956 [0.926, 1.074] 0.936 [0.905, 1.087] no misfit no misfit 0.40
PDSS.SR_7 0.915 [0.932, 1.085] 0.916 [0.921, 1.08] 0.017 0.005 0.22
Note:
MSQ values based on conditional calculations (n = 1678 complete cases).
Simulation based thresholds from 300 simulated datasets.
Code
RIgetfitPlot(simfit1, df_pre)

Code
ICCplot(as.data.frame(df_pre), 
        itemnumber = 1:4, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[1:4])

[1] "Please press Zoom on the Plots window to see the plot"
Code
### also suggested:
# library(RASCHplot) # install first with `devtools::install_github("ERRTG/RASCHplot")`
# CICCplot(PCM(df),
#          which.item = c(1:4),
#          lower.groups = c(0,7,14,21,28,35),
#          grid.items = TRUE)
Code
ICCplot(as.data.frame(df_pre), 
        itemnumber = 5:7, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[5:7])

[1] "Please press Zoom on the Plots window to see the plot"
Code
RIrestscore(df_pre)
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
PDSS.SR_1 0.63 0.58 0.05 0.011 * -0.10 0.62
PDSS.SR_2 0.66 0.59 0.07 0.000 *** -0.07 0.64
PDSS.SR_3 0.60 0.59 0.01 0.311 -0.96 -0.25
PDSS.SR_4 0.58 0.59 0.01 0.715 -0.24 0.47
PDSS.SR_5 0.47 0.60 0.13 0.000 *** 0.75 1.46
PDSS.SR_6 0.62 0.59 0.03 0.116 0.40 1.11
PDSS.SR_7 0.64 0.59 0.05 0.011 * 0.22 0.93
Code
RIbootRestscore(df_pre,cpu=n_cores,iterations = 500,samplesize=800) # samplesize=600
Item Item-restscore result % of iterations Conditional MSQ infit Relative average item location
PDSS.SR_1 no misfit 62.8 1.00 0.62
PDSS.SR_1 overfit 37.2 1.00 0.62
PDSS.SR_2 no misfit 16.6 0.87 0.64
PDSS.SR_2 overfit 83.4 0.87 0.64
PDSS.SR_3 no misfit 92.6 0.97 -0.25
PDSS.SR_3 overfit 7.4 0.97 -0.25
PDSS.SR_4 no misfit 99.4 1.05 0.47
PDSS.SR_5 underfit 99.8 1.29 1.46
PDSS.SR_6 no misfit 85.6 0.96 1.11
PDSS.SR_6 overfit 14.4 0.96 1.11
PDSS.SR_7 no misfit 63.2 0.91 0.93
PDSS.SR_7 overfit 36.8 0.91 0.93
Note:
Results based on 500 bootstrap iterations with a sample size of 800. Conditional mean-square infit based on complete responders only, n = 1678.
Code
RIpcmPCA(df_pre)
PCA of Rasch model residuals
Eigenvalues Proportion of variance
1.90 26.1%
1.28 21.3%
1.21 17.1%
1.07 14.7%
0.89 12.2%
Code
simcor1 <- RIgetResidCor(df_pre, iterations = 500, cpu = n_cores)
RIresidcorr(df_pre, cutoff = simcor1$p99)
PDSS.SR_1 PDSS.SR_2 PDSS.SR_3 PDSS.SR_4 PDSS.SR_5 PDSS.SR_6 PDSS.SR_7
PDSS.SR_1
PDSS.SR_2 0.37
PDSS.SR_3 -0.08 -0.05
PDSS.SR_4 -0.32 -0.28 -0.08
PDSS.SR_5 -0.24 -0.23 -0.22 -0.15
PDSS.SR_6 -0.15 -0.17 -0.2 -0.17 -0.14
PDSS.SR_7 -0.26 -0.26 -0.16 0.09 -0.2 -0.03
Note:
Relative cut-off value is -0.046, which is 0.094 above the average correlation (-0.14).
Correlations above the cut-off are highlighted in red text.
Code
RIpartgamLD(df_pre)
Item 1 Item 2 Partial gamma SE Lower CI Upper CI Adjusted p-value (BH)
PDSS.SR_2 PDSS.SR_1 0.737 0.026 0.686 0.787 0.000
PDSS.SR_1 PDSS.SR_2 0.724 0.027 0.670 0.778 0.000
PDSS.SR_7 PDSS.SR_4 0.34 0.038 0.266 0.415 0.000
PDSS.SR_4 PDSS.SR_7 0.305 0.038 0.230 0.380 0.000
PDSS.SR_7 PDSS.SR_6 0.212 0.039 0.135 0.288 0.000
PDSS.SR_6 PDSS.SR_7 0.187 0.040 0.108 0.265 0.000
PDSS.SR_2 PDSS.SR_3 0.164 0.043 0.080 0.248 0.006
PDSS.SR_3 PDSS.SR_4 0.141 0.041 0.062 0.221 0.021
Code
RIloadLoc(df_pre)

Code
mirt(df_pre, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Code
# for fewer items or a more magnified figure, use:
#RIitemCats(df)
Code
# increase fig-height above as needed, if you have many items
RItargeting(df_pre)

Code

Code
iarm::score_groups(as.data.frame(df_pre)) %>% 
  as.data.frame(nm = "score_group") %>% 
  dplyr::count(score_group)
  score_group   n
1           1 894
2           2 784
Code
dif_plots <- df_pre %>% 
  add_column(dif = iarm::score_groups(.)) %>% 
  split(.$dif) %>% # split the data using the DIF variable
  map(~ RItileplot(.x %>% dplyr::select(!dif)) + labs(title = .x$dif))
dif_plots[[1]] + dif_plots[[2]]

Code
clr_tests(df_pre, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr df pvalue  sig 
overall 115 27 7.4e-13  ***
Code
item_obsexp(PCM(df_pre))
Score group 1: 
          mean obs mean exp std.res sig
PDSS.SR_1  1.1466   1.1648  -0.9382    
PDSS.SR_2  1.3078   1.3388  -1.4004    
PDSS.SR_3  1.5896   1.5887   0.0406    
PDSS.SR_4  1.2255   1.2332  -0.3487    
PDSS.SR_5  0.6742   0.5928   3.8498 ++ 
PDSS.SR_6  0.7835   0.7779   0.2626    
PDSS.SR_7  0.9267   0.9450  -0.8361    

Score group 2: 
          mean obs mean exp std.res sig
PDSS.SR_1  2.107    2.087    0.882     
PDSS.SR_2  2.354    2.317    1.701     
PDSS.SR_3  2.696    2.685    0.466     
PDSS.SR_4  2.365    2.357    0.360     
PDSS.SR_5  1.801    1.893   -3.443  -- 
PDSS.SR_6  2.042    2.048   -0.238     
PDSS.SR_7  2.215    2.193    0.841     
Code
grouping_based_on_score = score_groups(as.data.frame(df_pre))
partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = grouping_based_on_score) 
       Item                     Var gamma  se pvalue padj.BH sig lower upper
1 PDSS.SR_1 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
2 PDSS.SR_2 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
3 PDSS.SR_3 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
4 PDSS.SR_4 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
5 PDSS.SR_5 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
6 PDSS.SR_6 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
7 PDSS.SR_7 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
Code
dif.sex <- dif %>% filter(Time=='PRE') %>% select(sex) %>% as.vector(.)
dif.sex <- dif.sex[[1]] #female = 1, male = 2 
RIdifTable(df_pre, dif.sex)
[1] "No statistically significant DIF found."
Code
partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.sex) 
       Item     Var   gamma     se pvalue padj.BH sig   lower   upper
1 PDSS.SR_1 dif.sex  0.0470 0.0522 0.3676  1.0000     -0.0553  0.1494
2 PDSS.SR_2 dif.sex  0.0564 0.0519 0.2777  1.0000     -0.0454  0.1582
3 PDSS.SR_3 dif.sex -0.0701 0.0489 0.1519  1.0000     -0.1661  0.0258
4 PDSS.SR_4 dif.sex -0.1046 0.0487 0.0316  0.2215     -0.2001 -0.0092
5 PDSS.SR_5 dif.sex  0.0693 0.0456 0.1280  0.8959     -0.0199  0.1586
6 PDSS.SR_6 dif.sex  0.0318 0.0483 0.5105  1.0000     -0.0629  0.1265
7 PDSS.SR_7 dif.sex -0.0061 0.0490 0.9015  1.0000     -0.1020  0.0899
Code
dif.age <- dif %>% filter(Time=='PRE') %>% select(age) %>% as.vector(.)
dif.age <- dif.age[[1]]
RIdifTable(df_pre, dif.age)

Item 2 3 Mean location StDev MaxDiff
PDSS.SR_1 -0.138 -0.001 -0.069 0.096 0.136
PDSS.SR_2 -0.108 0.070 -0.019 0.125 0.177
PDSS.SR_3 -1.122 -0.649 -0.885 0.335 0.474
PDSS.SR_4 -0.164 -0.445 -0.304 0.199 0.282
PDSS.SR_5 0.792 0.659 0.726 0.094 0.133
PDSS.SR_6 0.472 0.243 0.357 0.162 0.229
PDSS.SR_7 0.267 0.123 0.195 0.102 0.144
Code
partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.age) 
       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 PDSS.SR_1 dif.age -0.0639 0.0311 0.0402  0.2813      -0.1248 -0.0029
2 PDSS.SR_2 dif.age -0.0420 0.0305 0.1693  1.0000      -0.1018  0.0179
3 PDSS.SR_3 dif.age -0.1140 0.0278 0.0000  0.0003  *** -0.1686 -0.0595
4 PDSS.SR_4 dif.age  0.0919 0.0291 0.0016  0.0109    *  0.0350  0.1489
5 PDSS.SR_5 dif.age  0.0096 0.0277 0.7281  1.0000      -0.0447  0.0640
6 PDSS.SR_6 dif.age  0.0401 0.0277 0.1476  1.0000      -0.0142  0.0943
7 PDSS.SR_7 dif.age  0.0438 0.0287 0.1260  0.8817      -0.0123  0.1000
Code
dif.edu <- dif %>% filter(Time=='PRE') %>% select(education) %>% as.vector(.)
dif.edu <- dif.edu[[1]]
RIdifTable(df_pre, dif.edu)

Item 2 3 Mean location StDev MaxDiff
PDSS.SR_1 -0.078 -0.192 -0.135 0.081 0.114
PDSS.SR_2 -0.116 0.200 0.042 0.223 0.316
PDSS.SR_3 -0.908 -1.089 -0.999 0.128 0.181
PDSS.SR_4 -0.236 -0.304 -0.270 0.048 0.068
PDSS.SR_5 0.769 0.768 0.768 0.001 0.001
PDSS.SR_6 0.396 0.383 0.389 0.009 0.013
PDSS.SR_7 0.172 0.234 0.203 0.043 0.061
Code
partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.edu) 
       Item     Var   gamma     se pvalue padj.BH sig   lower   upper
1 PDSS.SR_1 dif.edu  0.0646 0.0470 0.1688  1.0000     -0.0274  0.1567
2 PDSS.SR_2 dif.edu  0.0282 0.0465 0.5443  1.0000     -0.0630  0.1195
3 PDSS.SR_3 dif.edu  0.0292 0.0435 0.5025  1.0000     -0.0561  0.1144
4 PDSS.SR_4 dif.edu -0.0733 0.0432 0.0900  0.6297     -0.1580  0.0114
5 PDSS.SR_5 dif.edu  0.1041 0.0410 0.0112  0.0783   .  0.0237  0.1844
6 PDSS.SR_6 dif.edu -0.0347 0.0433 0.4225  1.0000     -0.1196  0.0502
7 PDSS.SR_7 dif.edu -0.1511 0.0431 0.0005  0.0032  ** -0.2357 -0.0666
Code
dif.yr <- dif %>% filter(Time=='PRE') %>% select(TreatmentAccessStart) %>% as.vector(.)
dif.yr <- dif.yr[[1]]
RIdifTable(df_pre, dif.yr)
[1] "No statistically significant DIF found."
Code
partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.yr) 
       Item    Var   gamma     se pvalue padj.BH sig   lower  upper
1 PDSS.SR_1 dif.yr  0.0317 0.0329 0.3354  1.0000     -0.0328 0.0963
2 PDSS.SR_2 dif.yr -0.0079 0.0334 0.8120  1.0000     -0.0733 0.0574
3 PDSS.SR_3 dif.yr  0.0320 0.0296 0.2796  1.0000     -0.0260 0.0899
4 PDSS.SR_4 dif.yr -0.0377 0.0308 0.2206  1.0000     -0.0979 0.0226
5 PDSS.SR_5 dif.yr -0.0440 0.0277 0.1117  0.7819     -0.0983 0.0102
6 PDSS.SR_6 dif.yr  0.0475 0.0306 0.1206  0.8442     -0.0125 0.1075
7 PDSS.SR_7 dif.yr -0.0204 0.0301 0.4984  1.0000     -0.0795 0.0387

8.3 Analysis part 1 decision

From conditional item fit item 2 and 5 are flagged. However Item 2 have large residual correlations with 1 (possibly also seen at the contrast loadings), however item 5 is most different from the score groups. Removing item 5 due to the underfit.

9 Rasch analysis 2

Based on the above we remove item 5

Code
df_pre <- df %>% filter(Time=='PRE') %>% select(!c(Time,PDSS.SR_5)) #%>% sample_n(800)
RIlistItemsMargin(df_pre, fontsize = 12)
itemnr item
PDSS.SR_1 Number
PDSS.SR_2 Distress
PDSS.SR_3 Worry/Fear
PDSS.SR_4 Avoidance places
PDSS.SR_6 Impairment work/home
PDSS.SR_7 Impairment social
Code
RIitemfit(df_pre, cutoff = "Smith98")
Item InfitMSQ Location 1 ± 2 / √n
PDSS.SR_1 1.015 0.04 [0.951, 1.049]
PDSS.SR_2 0.897 0.04 [0.951, 1.049]
PDSS.SR_3 1.013 -0.88 [0.951, 1.049]
PDSS.SR_4 1.12 -0.12 [0.951, 1.049]
PDSS.SR_6 1.043 0.56 [0.951, 1.049]
PDSS.SR_7 0.961 0.36 [0.951, 1.049]
Note:
MSQ values based on conditional estimation (n = 1678 complete cases).
Code
simfit1 <- RIgetfit(df_pre, iterations = 300, cpu = n_cores)  

RIitemfit(df_pre, simfit1)
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Location
PDSS.SR_1 1.015 [0.93, 1.074] 0.992 [0.929, 1.07] no misfit no misfit 0.04
PDSS.SR_2 0.897 [0.922, 1.069] 0.896 [0.92, 1.071] 0.025 0.024 0.04
PDSS.SR_3 1.013 [0.931, 1.074] 1.012 [0.93, 1.079] no misfit no misfit -0.88
PDSS.SR_4 1.12 [0.93, 1.071] 1.132 [0.93, 1.071] 0.049 0.061 -0.12
PDSS.SR_6 1.043 [0.925, 1.096] 1.017 [0.918, 1.107] no misfit no misfit 0.56
PDSS.SR_7 0.961 [0.931, 1.079] 0.968 [0.919, 1.089] no misfit no misfit 0.36
Note:
MSQ values based on conditional calculations (n = 1678 complete cases).
Simulation based thresholds from 300 simulated datasets.
Code
RIgetfitPlot(simfit1, df_pre)

Code
ICCplot(as.data.frame(df_pre), 
        itemnumber = 1:4, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[1:4])

[1] "Please press Zoom on the Plots window to see the plot"
Code
### also suggested:
# library(RASCHplot) # install first with `devtools::install_github("ERRTG/RASCHplot")`
# CICCplot(PCM(df),
#          which.item = c(1:4),
#          lower.groups = c(0,7,14,21,28,35),
#          grid.items = TRUE)
Code
ICCplot(as.data.frame(df_pre), 
        itemnumber = 5:6, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[5:6])

[1] "Please press Zoom on the Plots window to see the plot"
Code
RIrestscore(df_pre)
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
PDSS.SR_1 0.64 0.61 0.03 0.189 0.04 0.65
PDSS.SR_2 0.67 0.62 0.05 0.003 ** 0.04 0.66
PDSS.SR_3 0.61 0.61 0.00 0.868 -0.88 -0.26
PDSS.SR_4 0.58 0.62 0.04 0.117 -0.12 0.49
PDSS.SR_6 0.61 0.62 0.01 0.697 0.56 1.17
PDSS.SR_7 0.64 0.62 0.02 0.282 0.36 0.97
Code
RIbootRestscore(df_pre,cpu=n_cores,iterations = 500,samplesize=800) # samplesize=600
Item Item-restscore result % of iterations Conditional MSQ infit Relative average item location
PDSS.SR_1 no misfit 94.8 1.02 0.65
PDSS.SR_1 overfit 5.2 1.02 0.65
PDSS.SR_2 no misfit 52.8 0.90 0.66
PDSS.SR_2 overfit 47.2 0.90 0.66
PDSS.SR_3 no misfit 99.4 1.01 -0.26
PDSS.SR_4 no misfit 91.0 1.12 0.49
PDSS.SR_4 underfit 9.0 1.12 0.49
PDSS.SR_6 no misfit 99.8 1.04 1.17
PDSS.SR_7 no misfit 96.4 0.96 0.97
Note:
Results based on 500 bootstrap iterations with a sample size of 800. Conditional mean-square infit based on complete responders only, n = 1678.
Code
RIpcmPCA(df_pre)
PCA of Rasch model residuals
Eigenvalues Proportion of variance
1.96 32.4%
1.30 22.1%
1.11 18.8%
0.93 15.8%
0.68 10.7%
Code
simcor1 <- RIgetResidCor(df_pre, iterations = 500, cpu = n_cores)
RIresidcorr(df_pre, cutoff = simcor1$p99)
PDSS.SR_1 PDSS.SR_2 PDSS.SR_3 PDSS.SR_4 PDSS.SR_6 PDSS.SR_7
PDSS.SR_1
PDSS.SR_2 0.33
PDSS.SR_3 -0.13 -0.11
PDSS.SR_4 -0.38 -0.34 -0.12
PDSS.SR_6 -0.2 -0.22 -0.25 -0.2
PDSS.SR_7 -0.33 -0.32 -0.22 0.05 -0.07
Note:
Relative cut-off value is -0.085, which is 0.083 above the average correlation (-0.167).
Correlations above the cut-off are highlighted in red text.
Code
RIpartgamLD(df_pre)
Item 1 Item 2 Partial gamma SE Lower CI Upper CI Adjusted p-value (BH)
PDSS.SR_2 PDSS.SR_1 0.737 0.026 0.686 0.788 0
PDSS.SR_1 PDSS.SR_2 0.707 0.028 0.652 0.763 0
PDSS.SR_7 PDSS.SR_4 0.352 0.037 0.279 0.426 0
PDSS.SR_4 PDSS.SR_7 0.303 0.038 0.228 0.379 0
PDSS.SR_7 PDSS.SR_6 0.235 0.039 0.158 0.312 0
PDSS.SR_6 PDSS.SR_7 0.179 0.041 0.100 0.259 0
Code
RIloadLoc(df_pre)

Code
mirt(df_pre, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Code
# for fewer items or a more magnified figure, use:
#RIitemCats(df)
Code
# increase fig-height above as needed, if you have many items
RItargeting(df_pre)

Code

Code
iarm::score_groups(as.data.frame(df_pre)) %>% 
  as.data.frame(nm = "score_group") %>% 
  dplyr::count(score_group)
  score_group   n
1           1 913
2           2 765
Code
dif_plots <- df_pre %>% 
  add_column(dif = iarm::score_groups(.)) %>% 
  split(.$dif) %>% # split the data using the DIF variable
  map(~ RItileplot(.x %>% dplyr::select(!dif)) + labs(title = .x$dif))
dif_plots[[1]] + dif_plots[[2]]

Code
clr_tests(df_pre, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr df pvalue  sig 
overall 60  23 3.8e-05  ***
Code
item_obsexp(PCM(df_pre))
Score group 1: 
          mean obs mean exp std.res sig
PDSS.SR_1  1.151    1.154   -0.171     
PDSS.SR_2  1.319    1.333   -0.658     
PDSS.SR_3  1.600    1.581    0.944     
PDSS.SR_4  1.236    1.224    0.597     
PDSS.SR_6  0.779    0.766    0.665     
PDSS.SR_7  0.914    0.934   -0.952     

Score group 2: 
          mean obs mean exp std.res sig
PDSS.SR_1  2.123    2.119    0.164     
PDSS.SR_2  2.364    2.347    0.828     
PDSS.SR_3  2.709    2.723   -0.639     
PDSS.SR_4  2.379    2.394   -0.629     
PDSS.SR_6  2.076    2.092   -0.618     
PDSS.SR_7  2.260    2.236    0.961     
Code
grouping_based_on_score = score_groups(as.data.frame(df_pre))
partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = grouping_based_on_score) 
       Item                     Var gamma  se pvalue padj.BH sig lower upper
1 PDSS.SR_1 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
2 PDSS.SR_2 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
3 PDSS.SR_3 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
4 PDSS.SR_4 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
5 PDSS.SR_6 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
6 PDSS.SR_7 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
Code
dif.sex <- dif %>% filter(Time=='PRE') %>% select(sex) %>% as.vector(.)
dif.sex <- dif.sex[[1]] #female = 1, male = 2 
RIdifTable(df_pre, dif.sex)
[1] "No statistically significant DIF found."
Code
partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.sex) 
       Item     Var   gamma     se pvalue padj.BH sig   lower  upper
1 PDSS.SR_1 dif.sex  0.0540 0.0534 0.3118  1.0000     -0.0507 0.1587
2 PDSS.SR_2 dif.sex  0.0513 0.0526 0.3290  1.0000     -0.0517 0.1543
3 PDSS.SR_3 dif.sex -0.0725 0.0496 0.1437  0.8625     -0.1698 0.0247
4 PDSS.SR_4 dif.sex -0.0911 0.0489 0.0627  0.3764     -0.1869 0.0048
5 PDSS.SR_6 dif.sex  0.0572 0.0484 0.2374  1.0000     -0.0377 0.1522
6 PDSS.SR_7 dif.sex  0.0195 0.0493 0.6924  1.0000     -0.0772 0.1162
Code
dif.age <- dif %>% filter(Time=='PRE') %>% select(age) %>% as.vector(.)
dif.age <- dif.age[[1]]
RIdifTable(df_pre, dif.age)

Item 2 3 Mean location StDev MaxDiff
PDSS.SR_1 0.004 0.124 0.064 0.085 0.120
PDSS.SR_2 0.008 0.175 0.092 0.118 0.167
PDSS.SR_3 -1.030 -0.574 -0.802 0.322 0.456
PDSS.SR_4 -0.032 -0.362 -0.197 0.234 0.330
PDSS.SR_6 0.634 0.390 0.512 0.172 0.244
PDSS.SR_7 0.416 0.247 0.332 0.119 0.168
Code
partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.age) 
       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 PDSS.SR_1 dif.age -0.0643 0.0316 0.0414  0.2486      -0.1262 -0.0025
2 PDSS.SR_2 dif.age -0.0418 0.0310 0.1780  1.0000      -0.1026  0.0190
3 PDSS.SR_3 dif.age -0.1211 0.0280 0.0000  0.0001  *** -0.1760 -0.0662
4 PDSS.SR_4 dif.age  0.1013 0.0292 0.0005  0.0031   **  0.0441  0.1584
5 PDSS.SR_6 dif.age  0.0354 0.0278 0.2023  1.0000      -0.0190  0.0899
6 PDSS.SR_7 dif.age  0.0485 0.0287 0.0907  0.5439      -0.0077  0.1047
Code
dif.edu <- dif %>% filter(Time=='PRE') %>% select(education) %>% as.vector(.)
dif.edu <- dif.edu[[1]]
RIdifTable(df_pre, dif.edu)
[1] "No statistically significant DIF found."
Code
partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.edu) 
       Item     Var   gamma     se pvalue padj.BH sig   lower   upper
1 PDSS.SR_1 dif.edu  0.1017 0.0473 0.0316  0.1894      0.0090  0.1945
2 PDSS.SR_2 dif.edu  0.0683 0.0471 0.1467  0.8800     -0.0239  0.1606
3 PDSS.SR_3 dif.edu  0.0585 0.0438 0.1815  1.0000     -0.0273  0.1442
4 PDSS.SR_4 dif.edu -0.0607 0.0436 0.1637  0.9821     -0.1462  0.0247
5 PDSS.SR_6 dif.edu -0.0150 0.0432 0.7288  1.0000     -0.0996  0.0696
6 PDSS.SR_7 dif.edu -0.1315 0.0438 0.0027  0.0162   * -0.2174 -0.0456
Code
dif.yr <- dif %>% filter(Time=='PRE') %>% select(TreatmentAccessStart) %>% as.vector(.)
dif.yr <- dif.yr[[1]]
RIdifTable(df_pre, dif.yr)
[1] "No statistically significant DIF found."
Code
partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.yr) 
       Item    Var   gamma     se pvalue padj.BH sig   lower  upper
1 PDSS.SR_1 dif.yr  0.0317 0.0330 0.3369  1.0000     -0.0330 0.0963
2 PDSS.SR_2 dif.yr -0.0090 0.0337 0.7901  1.0000     -0.0750 0.0571
3 PDSS.SR_3 dif.yr  0.0307 0.0301 0.3069  1.0000     -0.0282 0.0896
4 PDSS.SR_4 dif.yr -0.0525 0.0306 0.0864  0.5184     -0.1126 0.0075
5 PDSS.SR_6 dif.yr  0.0427 0.0302 0.1569  0.9414     -0.0164 0.1019
6 PDSS.SR_7 dif.yr -0.0357 0.0305 0.2411  1.0000     -0.0954 0.0240

9.1 Analysis part 2 decision

Item 2 displays the most misfit and most residual correlations (with 1 too).

10 Rasch analysis 3

Based on the above we remove item 2

Code
df_pre <- df %>% filter(Time=='PRE') %>% select(!c(Time,PDSS.SR_5,PDSS.SR_2)) #%>% sample_n(800)
RIlistItemsMargin(df_pre, fontsize = 12)
itemnr item
PDSS.SR_1 Number
PDSS.SR_3 Worry/Fear
PDSS.SR_4 Avoidance places
PDSS.SR_6 Impairment work/home
PDSS.SR_7 Impairment social
Code
RIitemfit(df_pre, cutoff = "Smith98")
Item InfitMSQ Location 1 ± 2 / √n
PDSS.SR_1 1.151 0.04 [0.951, 1.049]
PDSS.SR_3 1.002 -0.84 [0.951, 1.049]
PDSS.SR_4 1.02 -0.11 [0.951, 1.049]
PDSS.SR_6 0.989 0.55 [0.951, 1.049]
PDSS.SR_7 0.869 0.36 [0.951, 1.049]
Note:
MSQ values based on conditional estimation (n = 1678 complete cases).
Code
simfit1 <- RIgetfit(df_pre, iterations = 300, cpu = n_cores)  

RIitemfit(df_pre, simfit1)
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Location
PDSS.SR_1 1.151 [0.927, 1.088] 1.141 [0.928, 1.086] 0.063 0.055 0.04
PDSS.SR_3 1.002 [0.928, 1.069] 0.996 [0.931, 1.068] no misfit no misfit -0.84
PDSS.SR_4 1.02 [0.935, 1.063] 1.023 [0.933, 1.075] no misfit no misfit -0.11
PDSS.SR_6 0.989 [0.938, 1.077] 0.969 [0.931, 1.101] no misfit no misfit 0.55
PDSS.SR_7 0.869 [0.924, 1.054] 0.864 [0.922, 1.065] 0.055 0.058 0.36
Note:
MSQ values based on conditional calculations (n = 1678 complete cases).
Simulation based thresholds from 300 simulated datasets.
Code
RIgetfitPlot(simfit1, df_pre)

Code
ICCplot(as.data.frame(df_pre), 
        itemnumber = 1:4, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[c(1,3,4,6)])

[1] "Please press Zoom on the Plots window to see the plot"
Code
### also suggested:
# library(RASCHplot) # install first with `devtools::install_github("ERRTG/RASCHplot")`
# CICCplot(PCM(df),
#          which.item = c(1:4),
#          lower.groups = c(0,7,14,21,28,35),
#          grid.items = TRUE)
Code
ICCplot(as.data.frame(df_pre), 
        itemnumber = 5, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[7])

Code
RIrestscore(df_pre)
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
PDSS.SR_1 0.55 0.59 0.04 0.068 . 0.04 0.63
PDSS.SR_3 0.59 0.60 0.01 0.962 -0.84 -0.26
PDSS.SR_4 0.60 0.60 0.00 0.962 -0.11 0.47
PDSS.SR_6 0.61 0.60 0.01 0.883 0.55 1.14
PDSS.SR_7 0.67 0.60 0.07 0.000 *** 0.36 0.94
Code
RIbootRestscore(df_pre,cpu=n_cores,iterations = 500,samplesize=800) # samplesize=600
Item Item-restscore result % of iterations Conditional MSQ infit Relative average item location
PDSS.SR_1 no misfit 88.4 1.15 0.63
PDSS.SR_1 underfit 11.6 1.15 0.63
PDSS.SR_3 no misfit 99.8 1.00 -0.26
PDSS.SR_4 no misfit 99.6 1.02 0.47
PDSS.SR_6 no misfit 97.8 0.99 1.14
PDSS.SR_7 no misfit 34.0 0.87 0.94
PDSS.SR_7 overfit 66.0 0.87 0.94
Note:
Results based on 500 bootstrap iterations with a sample size of 800. Conditional mean-square infit based on complete responders only, n = 1678.
Code
RIpcmPCA(df_pre)
PCA of Rasch model residuals
Eigenvalues Proportion of variance
1.55 32.2%
1.37 26.9%
1.06 21.7%
1.01 19%
0.01 0.1%
Code
simcor1 <- RIgetResidCor(df_pre, iterations = 500, cpu = n_cores)
RIresidcorr(df_pre, cutoff = simcor1$p99)
PDSS.SR_1 PDSS.SR_3 PDSS.SR_4 PDSS.SR_6 PDSS.SR_7
PDSS.SR_1
PDSS.SR_3 -0.07
PDSS.SR_4 -0.35 -0.17
PDSS.SR_6 -0.15 -0.28 -0.28
PDSS.SR_7 -0.29 -0.27 -0.03 -0.13
Note:
Relative cut-off value is -0.116, which is 0.085 above the average correlation (-0.201).
Correlations above the cut-off are highlighted in red text.
Code
RIpartgamLD(df_pre)
Item 1 Item 2 Partial gamma SE Lower CI Upper CI Adjusted p-value (BH)
PDSS.SR_7 PDSS.SR_4 0.314 0.039 0.237 0.391 0.000
PDSS.SR_4 PDSS.SR_7 0.224 0.041 0.144 0.304 0.000
PDSS.SR_3 PDSS.SR_1 0.203 0.041 0.123 0.284 0.000
PDSS.SR_7 PDSS.SR_6 0.197 0.041 0.117 0.277 0.000
PDSS.SR_1 PDSS.SR_3 0.16 0.041 0.079 0.241 0.002
PDSS.SR_6 PDSS.SR_1 0.138 0.042 0.055 0.220 0.021
Code
RIloadLoc(df_pre)

Code
mirt(df_pre, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Code
# for fewer items or a more magnified figure, use:
#RIitemCats(df)
Code
# increase fig-height above as needed, if you have many items
RItargeting(df_pre)

Code

Code
iarm::score_groups(as.data.frame(df_pre)) %>% 
  as.data.frame(nm = "score_group") %>% 
  dplyr::count(score_group)
  score_group   n
1           1 896
2           2 782
Code
dif_plots <- df_pre %>% 
  add_column(dif = iarm::score_groups(.)) %>% 
  split(.$dif) %>% # split the data using the DIF variable
  map(~ RItileplot(.x %>% dplyr::select(!dif)) + labs(title = .x$dif))
dif_plots[[1]] + dif_plots[[2]]

Code
clr_tests(df_pre, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr  df pvalue  sig 
overall 56.7 19 1.3e-05  ***
Code
item_obsexp(PCM(df_pre))
Score group 1: 
          mean obs mean exp std.res sig
PDSS.SR_1  1.170    1.146    1.283     
PDSS.SR_3  1.586    1.570    0.768     
PDSS.SR_4  1.207    1.212   -0.211     
PDSS.SR_6  0.761    0.752    0.478     
PDSS.SR_7  0.882    0.919   -1.829     

Score group 2: 
          mean obs mean exp std.res sig
PDSS.SR_1  2.082    2.109   -1.198     
PDSS.SR_3  2.704    2.714   -0.471     
PDSS.SR_4  2.388    2.383    0.215     
PDSS.SR_6  2.069    2.080   -0.419     
PDSS.SR_7  2.268    2.225    1.761     
Code
grouping_based_on_score = score_groups(as.data.frame(df_pre))
partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = grouping_based_on_score) 
       Item                     Var gamma  se pvalue padj.BH sig lower upper
1 PDSS.SR_1 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
2 PDSS.SR_3 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
3 PDSS.SR_4 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
4 PDSS.SR_6 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
5 PDSS.SR_7 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
Code
dif.sex <- dif %>% filter(Time=='PRE') %>% select(sex) %>% as.vector(.)
dif.sex <- dif.sex[[1]] #female = 1, male = 2 
RIdifTable(df_pre, dif.sex)
[1] "No statistically significant DIF found."
Code
partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.sex) 
       Item     Var   gamma     se pvalue padj.BH sig   lower  upper
1 PDSS.SR_1 dif.sex  0.0469 0.0506 0.3544  1.0000     -0.0523 0.1461
2 PDSS.SR_3 dif.sex -0.0569 0.0494 0.2490  1.0000     -0.1537 0.0398
3 PDSS.SR_4 dif.sex -0.0863 0.0496 0.0818  0.4089     -0.1835 0.0109
4 PDSS.SR_6 dif.sex  0.0618 0.0489 0.2057  1.0000     -0.0339 0.1576
5 PDSS.SR_7 dif.sex  0.0364 0.0505 0.4707  1.0000     -0.0625 0.1353
Code
dif.age <- dif %>% filter(Time=='PRE') %>% select(age) %>% as.vector(.)
dif.age <- dif.age[[1]]
RIdifTable(df_pre, dif.age)

Item 2 3 Mean location StDev MaxDiff
PDSS.SR_1 0.000 0.149 0.075 0.105 0.149
PDSS.SR_3 -1.003 -0.525 -0.764 0.338 0.478
PDSS.SR_4 -0.028 -0.313 -0.171 0.201 0.285
PDSS.SR_6 0.621 0.410 0.516 0.149 0.211
PDSS.SR_7 0.409 0.278 0.344 0.093 0.131
Code
partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.age) 
       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 PDSS.SR_1 dif.age -0.0710 0.0302 0.0187  0.0933    . -0.1301 -0.0118
2 PDSS.SR_3 dif.age -0.1377 0.0277 0.0000  0.0000  *** -0.1920 -0.0834
3 PDSS.SR_4 dif.age  0.0991 0.0294 0.0008  0.0038   **  0.0414  0.1568
4 PDSS.SR_6 dif.age  0.0393 0.0284 0.1666  0.8332      -0.0164  0.0950
5 PDSS.SR_7 dif.age  0.0446 0.0294 0.1297  0.6483      -0.0131  0.1023
Code
dif.edu <- dif %>% filter(Time=='PRE') %>% select(education) %>% as.vector(.)
dif.edu <- dif.edu[[1]]
RIdifTable(df_pre, dif.edu)
[1] "No statistically significant DIF found."
Code
partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.edu) 
       Item     Var   gamma     se pvalue padj.BH sig   lower   upper
1 PDSS.SR_1 dif.edu  0.0824 0.0457 0.0711  0.3557     -0.0071  0.1719
2 PDSS.SR_3 dif.edu  0.0803 0.0438 0.0669  0.3344     -0.0056  0.1661
3 PDSS.SR_4 dif.edu -0.0374 0.0439 0.3942  1.0000     -0.1235  0.0486
4 PDSS.SR_6 dif.edu -0.0055 0.0441 0.9007  1.0000     -0.0920  0.0810
5 PDSS.SR_7 dif.edu -0.1295 0.0450 0.0040  0.0198   * -0.2176 -0.0414
Code
dif.yr <- dif %>% filter(Time=='PRE') %>% select(TreatmentAccessStart) %>% as.vector(.)
dif.yr <- dif.yr[[1]]
RIdifTable(df_pre, dif.yr)
[1] "No statistically significant DIF found."
Code
partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.yr) 
       Item    Var   gamma     se pvalue padj.BH sig   lower  upper
1 PDSS.SR_1 dif.yr  0.0310 0.0319 0.3319  1.0000     -0.0316 0.0935
2 PDSS.SR_3 dif.yr  0.0422 0.0297 0.1551  0.7754     -0.0160 0.1004
3 PDSS.SR_4 dif.yr -0.0599 0.0313 0.0552  0.2760     -0.1212 0.0013
4 PDSS.SR_6 dif.yr  0.0414 0.0305 0.1739  0.8695     -0.0183 0.1011
5 PDSS.SR_7 dif.yr -0.0573 0.0308 0.0626  0.3129     -0.1176 0.0030

10.1 Analysis part 3 decision

The only really underperforming item is item 7. With the overfit indicated by the boostrap and the association / local dependence between item 4 and 7 as indicated by the partial gamma.

11 Rasch analysis 4

Based on the above we remove item 7

Code
df_pre <- df %>% filter(Time=='PRE') %>% select(!c(Time,PDSS.SR_5,PDSS.SR_2,PDSS.SR_7)) #%>% sample_n(800)
RIlistItemsMargin(df_pre, fontsize = 12)
itemnr item
PDSS.SR_1 Number
PDSS.SR_3 Worry/Fear
PDSS.SR_4 Avoidance places
PDSS.SR_6 Impairment work/home
Code
RIitemfit(df_pre, cutoff = "Smith98")
Item InfitMSQ Location 1 ± 2 / √n
PDSS.SR_1 1.053 0.12 [0.951, 1.049]
PDSS.SR_3 0.92 -0.72 [0.951, 1.049]
PDSS.SR_4 1.07 -0.02 [0.951, 1.049]
PDSS.SR_6 0.993 0.62 [0.951, 1.049]
Note:
MSQ values based on conditional estimation (n = 1678 complete cases).
Code
simfit1 <- RIgetfit(df_pre, iterations = 300, cpu = n_cores)  

RIitemfit(df_pre, simfit1)
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Location
PDSS.SR_1 1.053 [0.916, 1.076] 1.043 [0.921, 1.077] no misfit no misfit 0.12
PDSS.SR_3 0.92 [0.931, 1.072] 0.918 [0.93, 1.077] 0.011 0.012 -0.72
PDSS.SR_4 1.07 [0.929, 1.074] 1.077 [0.927, 1.08] no misfit no misfit -0.02
PDSS.SR_6 0.993 [0.934, 1.064] 0.974 [0.915, 1.084] no misfit no misfit 0.62
Note:
MSQ values based on conditional calculations (n = 1678 complete cases).
Simulation based thresholds from 300 simulated datasets.
Code
RIgetfitPlot(simfit1, df_pre)

Code
ICCplot(as.data.frame(df_pre), 
        itemnumber = 1:4, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[c(1,3,4,6)])

[1] "Please press Zoom on the Plots window to see the plot"
Code
### also suggested:
# library(RASCHplot) # install first with `devtools::install_github("ERRTG/RASCHplot")`
# CICCplot(PCM(df),
#          which.item = c(1:4),
#          lower.groups = c(0,7,14,21,28,35),
#          grid.items = TRUE)
Code
RIrestscore(df_pre)
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
PDSS.SR_1 0.56 0.56 0.00 0.748 0.12 0.59
PDSS.SR_3 0.60 0.56 0.04 0.049 * -0.72 -0.25
PDSS.SR_4 0.54 0.56 0.02 0.462 -0.02 0.45
PDSS.SR_6 0.57 0.57 0.00 0.748 0.62 1.09
Code
RIbootRestscore(df_pre,cpu=n_cores,iterations = 500,samplesize=800) # samplesize=600
Item Item-restscore result % of iterations Conditional MSQ infit Relative average item location
PDSS.SR_1 no misfit 99.8 1.05 0.59
PDSS.SR_3 no misfit 83.8 0.92 -0.25
PDSS.SR_3 overfit 16.2 0.92 -0.25
PDSS.SR_4 no misfit 99.8 1.07 0.45
PDSS.SR_6 no misfit 100.0 0.99 1.09
Note:
Results based on 500 bootstrap iterations with a sample size of 800. Conditional mean-square infit based on complete responders only, n = 1678.
Code
RIpcmPCA(df_pre)
PCA of Rasch model residuals
Eigenvalues Proportion of variance
1.47 38.8%
1.41 33.4%
1.11 27.6%
0.01 0.2%
Code
simcor1 <- RIgetResidCor(df_pre, iterations = 500, cpu = n_cores)
RIresidcorr(df_pre, cutoff = simcor1$p99)
PDSS.SR_1 PDSS.SR_3 PDSS.SR_4 PDSS.SR_6
PDSS.SR_1
PDSS.SR_3 -0.14
PDSS.SR_4 -0.37 -0.18
PDSS.SR_6 -0.19 -0.32 -0.25
Note:
Relative cut-off value is -0.162, which is 0.081 above the average correlation (-0.242).
Correlations above the cut-off are highlighted in red text.
Code
RIpartgamLD(df_pre)
Item 1 Item 2 Partial gamma SE Lower CI Upper CI Adjusted p-value (BH)
PDSS.SR_3 PDSS.SR_4 0.175 0.042 0.093 0.257 0.000
PDSS.SR_3 PDSS.SR_1 0.152 0.044 0.065 0.238 0.007
PDSS.SR_6 PDSS.SR_1 0.13 0.043 0.045 0.215 0.031
PDSS.SR_1 PDSS.SR_6 0.122 0.042 0.039 0.205 0.049
Code
RIloadLoc(df_pre)

Code
mirt(df_pre, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Code
# for fewer items or a more magnified figure, use:
#RIitemCats(df)
Code
# increase fig-height above as needed, if you have many items
RItargeting(df_pre)

Code

Code
iarm::score_groups(as.data.frame(df_pre)) %>% 
  as.data.frame(nm = "score_group") %>% 
  dplyr::count(score_group)
  score_group   n
1           1 790
2           2 888
Code
dif_plots <- df_pre %>% 
  add_column(dif = iarm::score_groups(.)) %>% 
  split(.$dif) %>% # split the data using the DIF variable
  map(~ RItileplot(.x %>% dplyr::select(!dif)) + labs(title = .x$dif))
dif_plots[[1]] + dif_plots[[2]]

Code
clr_tests(df_pre, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr  df pvalue  sig 
overall 38.5 15 0.00076  ***
Code
item_obsexp(PCM(df_pre))
Score group 1: 
          mean obs mean exp std.res sig
PDSS.SR_1  1.069    1.089   -1.002     
PDSS.SR_3  1.485    1.493   -0.340     
PDSS.SR_4  1.150    1.131    0.864     
PDSS.SR_6  0.684    0.672    0.564     

Score group 2: 
          mean obs mean exp std.res sig
PDSS.SR_1  2.061    2.044    0.832     
PDSS.SR_3  2.659    2.650    0.441     
PDSS.SR_4  2.297    2.314   -0.764     
PDSS.SR_6  1.980    1.990   -0.437     
Code
grouping_based_on_score = score_groups(as.data.frame(df_pre))
partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = grouping_based_on_score) 
       Item                     Var gamma  se pvalue padj.BH sig lower upper
1 PDSS.SR_1 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
2 PDSS.SR_3 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
3 PDSS.SR_4 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
4 PDSS.SR_6 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
Code
dif.sex <- dif %>% filter(Time=='PRE') %>% select(sex) %>% as.vector(.)
dif.sex <- dif.sex[[1]] #female = 1, male = 2 
RIdifTable(df_pre, dif.sex)
[1] "No statistically significant DIF found."
Code
partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.sex) 
       Item     Var   gamma     se pvalue padj.BH sig   lower  upper
1 PDSS.SR_1 dif.sex  0.0694 0.0522 0.1841  0.7364     -0.0330 0.1718
2 PDSS.SR_3 dif.sex -0.0514 0.0511 0.3141  1.0000     -0.1516 0.0487
3 PDSS.SR_4 dif.sex -0.0816 0.0492 0.0969  0.3877     -0.1780 0.0147
4 PDSS.SR_6 dif.sex  0.0723 0.0489 0.1398  0.5592     -0.0237 0.1682
Code
dif.age <- dif %>% filter(Time=='PRE') %>% select(age) %>% as.vector(.)
dif.age <- dif.age[[1]]
RIdifTable(df_pre, dif.age)

Item 2 3 Mean location StDev MaxDiff
PDSS.SR_1 0.093 0.201 0.147 0.077 0.109
PDSS.SR_3 -0.867 -0.430 -0.648 0.309 0.437
PDSS.SR_4 0.073 -0.227 -0.077 0.212 0.299
PDSS.SR_6 0.701 0.455 0.578 0.174 0.246
Code
partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.age) 
       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 PDSS.SR_1 dif.age -0.0643 0.0312 0.0391  0.1563      -0.1253 -0.0032
2 PDSS.SR_3 dif.age -0.1356 0.0291 0.0000  0.0000  *** -0.1927 -0.0786
3 PDSS.SR_4 dif.age  0.1023 0.0287 0.0004  0.0014   **  0.0462  0.1585
4 PDSS.SR_6 dif.age  0.0487 0.0286 0.0885  0.3539      -0.0073  0.1047
Code
dif.edu <- dif %>% filter(Time=='PRE') %>% select(education) %>% as.vector(.)
dif.edu <- dif.edu[[1]]
RIdifTable(df_pre, dif.edu)
[1] "No statistically significant DIF found."
Code
partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.edu) 
       Item     Var   gamma     se pvalue padj.BH sig   lower  upper
1 PDSS.SR_1 dif.edu  0.0600 0.0469 0.2011  0.8042     -0.0320 0.1519
2 PDSS.SR_3 dif.edu  0.0518 0.0452 0.2511  1.0000     -0.0367 0.1403
3 PDSS.SR_4 dif.edu -0.0803 0.0434 0.0646  0.2584     -0.1654 0.0049
4 PDSS.SR_6 dif.edu -0.0279 0.0440 0.5259  1.0000     -0.1141 0.0583
Code
dif.yr <- dif %>% filter(Time=='PRE') %>% select(TreatmentAccessStart) %>% as.vector(.)
dif.yr <- dif.yr[[1]]
RIdifTable(df_pre, dif.yr)
[1] "No statistically significant DIF found."
Code
partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.yr) 
       Item    Var   gamma     se pvalue padj.BH sig   lower   upper
1 PDSS.SR_1 dif.yr  0.0158 0.0325 0.6261  1.0000     -0.0478  0.0794
2 PDSS.SR_3 dif.yr  0.0201 0.0311 0.5185  1.0000     -0.0409  0.0811
3 PDSS.SR_4 dif.yr -0.0661 0.0308 0.0321  0.1285     -0.1265 -0.0056
4 PDSS.SR_6 dif.yr  0.0277 0.0309 0.3688  1.0000     -0.0328  0.0883

11.1 Analysis part 4 decision

While item 3 have some small overfit in the boostrap (~15%) it is not so severe. As is the small DIF indications on age.

12 Analysis conclusion

Resulting PDSS-SR reworked scale thus has item 1, 3, 4 and 6.

12.1 Resulting items

Code
df_pre <- df %>% filter(Time=='PRE') %>% select(!c(Time,PDSS.SR_5,PDSS.SR_2,PDSS.SR_7)) #%>% sample_n(800)

what_scale <- gsub('\\.','',items_to_use)
Code
RIlistItemsMargin(df_pre, fontsize = 12)
itemnr item
PDSS.SR_1 Number
PDSS.SR_3 Worry/Fear
PDSS.SR_4 Avoidance places
PDSS.SR_6 Impairment work/home

12.2 Reliability

Code
RItif(df_pre)

Code
RItif(df_pre,samplePSI=TRUE)

12.3 Person parameter

Code
RIpfit(df_pre)

12.4 Misfit check

Code
pfit_u3poly <- PerFit::U3poly(matrix = df_pre, 
                      Ncat = 5, # make sure to input number of response categories, not thresholds
                      IRT.PModel = "PCM")
misfits <- PerFit::flagged.resp(pfit_u3poly) %>% 
  pluck("Scores") %>% 
  as.data.frame() %>% 
  pull(FlaggedID)

misfits2 <- RIpfit(df_pre, output = "rowid")
Code
RIitemfit(df_pre, simcut = simfit1)
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Location
PDSS.SR_1 1.053 [0.916, 1.076] 1.043 [0.921, 1.077] no misfit no misfit 0.12
PDSS.SR_3 0.92 [0.931, 1.072] 0.918 [0.93, 1.077] 0.011 0.012 -0.72
PDSS.SR_4 1.07 [0.929, 1.074] 1.077 [0.927, 1.08] no misfit no misfit -0.02
PDSS.SR_6 0.993 [0.934, 1.064] 0.974 [0.915, 1.084] no misfit no misfit 0.62
Note:
MSQ values based on conditional calculations (n = 1678 complete cases).
Simulation based thresholds from 300 simulated datasets.
Code
RIitemfit(df_pre[-misfits,], simcut = simfit1)
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Location
PDSS.SR_1 1.006 [0.916, 1.076] 0.991 [0.921, 1.077] no misfit no misfit 0.97
PDSS.SR_3 0.98 [0.931, 1.072] 0.937 [0.93, 1.077] no misfit no misfit -3.50
PDSS.SR_4 1.017 [0.929, 1.074] 1.012 [0.927, 1.08] no misfit no misfit 0.79
PDSS.SR_6 1.015 [0.934, 1.064] 0.977 [0.915, 1.084] no misfit no misfit 1.74
Note:
MSQ values based on conditional calculations (n = 1495 complete cases).
Simulation based thresholds from 300 simulated datasets.
Code
RItif(df_pre[-misfits,])

Code
RIitemfit(df_pre[-misfits2,], simcut = simfit1)
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Location
PDSS.SR_1 1.073 [0.916, 1.076] 1.065 [0.921, 1.077] no misfit no misfit 0.15
PDSS.SR_3 0.942 [0.931, 1.072] 0.935 [0.93, 1.077] no misfit no misfit -0.79
PDSS.SR_4 1.016 [0.929, 1.074] 1.023 [0.927, 1.08] no misfit no misfit -0.02
PDSS.SR_6 0.998 [0.934, 1.064] 0.978 [0.915, 1.084] no misfit no misfit 0.66
Note:
MSQ values based on conditional calculations (n = 1559 complete cases).
Simulation based thresholds from 300 simulated datasets.
Code
RItif(df_pre[-misfits2,])

12.5 Item parameters

Code
RIitemparams(df_pre)
Threshold 1 Threshold 2 Threshold 3 Threshold 4 Item location
PDSS.SR_1 -3.17 -0.54 1.88 2.30 0.12
PDSS.SR_3 -4.02 -1.40 0.11 2.42 -0.72
PDSS.SR_4 -2.70 -0.93 0.97 2.59 -0.02
PDSS.SR_6 -1.69 0.04 1.05 3.07 0.62
Note:
Item location is the average of the thresholds for each item.
Code
item_param_matrix <- RIitemparams(df_pre,output='dataframe') %>% as.matrix(.)
write.csv(item_param_matrix,file=paste0("./results/item_params_",what_scale,'.csv'),row.names=TRUE)

12.6 Ordinal sum to interval score

Code
RIscoreSE(df_pre)
Ordinal sum score Logit score Logit std.error
0 -5.419 0.775
1 -3.960 0.968
2 -3.072 0.901
3 -2.367 0.818
4 -1.770 0.760
5 -1.243 0.721
6 -0.759 0.696
7 -0.300 0.679
8 0.146 0.669
9 0.587 0.663
10 1.025 0.662
11 1.460 0.668
12 1.892 0.689
13 2.340 0.731
14 2.849 0.793
15 3.523 0.842
16 4.787 0.693
Code
RIscoreSE(df_pre,output='figure')

Code
sum_to_latent <- RIscoreSE(df_pre,output = 'dataframe')
write.csv(sum_to_latent,file=paste0('./results/ordinal_sum_to_latent_',what_scale,'.csv'),row.names=FALSE)

12.7 Thetas

The ordinal sum to interval score contains the information to transform into the below thetas, but we plot them here for convinience (based on 4 items)

Code
est_thetas2 <- RIestThetas(df_pre, method = "WLE")
hist(est_thetas2$WLE, 
     col = "#009ca6", 
     main = "Histogram of person locations (thetas)", 
     breaks = 10)

13 Software used

Code
pkgs <- grateful::cite_packages(cite.tidyverse = TRUE, 
                      output = "table",
                      bib.file = "grateful-refs.bib",
                      include.RStudio = FALSE,
                      out.dir = getwd())
formattable(pkgs, table.attr = 'class=\"table table-striped\" style="font-size: 15px; font-family: Lato; width: 80%"')
Package Version Citation
base 4.4.1 @base
car 3.1.3 @car
doParallel 1.0.17 @doParallel
easyRasch 0.3.3 @easyRasch
eRm 1.0.6 @eRm2007a; @eRm2007b; @eRm2009c; @eRm2013d; @eRm2015e; @eRm2019f
formattable 0.2.1 @formattable
furrr 0.3.1 @furrr
ggrepel 0.9.6 @ggrepel
glue 1.8.0 @glue
gridExtra 2.3 @gridExtra
iarm 0.4.3 @iarm
janitor 2.2.0 @janitor
kableExtra 1.4.0 @kableExtra
knitr 1.49 @knitr2014; @knitr2015; @knitr2024
matrixStats 1.4.1 @matrixStats
mirt 1.43 @mirt
mokken 3.1.2 @mokken2007; @mokken2012
patchwork 1.3.0 @patchwork
PerFit 1.4.6 @PerFit
psych 2.4.6.26 @psych
psychotools 0.7.4 @psychotools2021; @psychotools2022; @psychotools2024
psychotree 0.16.1 @psychotree2010e; @psychotree2011a; @psychotree2015b; @psychotree2018c; @psychotree2018d
reshape 0.8.9 @reshape
rmarkdown 2.29 @rmarkdown2018; @rmarkdown2020; @rmarkdown2024
tidyverse 2.0.0 @tidyverse

14 Session info

Code
R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Europe/Stockholm
tzcode source: system (glibc)

attached base packages:
 [1] parallel  grid      stats4    stats     graphics  grDevices utils    
 [8] datasets  methods   base     

other attached packages:
 [1] gridExtra_2.3     knitr_1.49        car_3.1-3         carData_3.0-5    
 [5] grateful_0.2.10   easyRasch_0.3.3   doParallel_1.0.17 iterators_1.0.14 
 [9] furrr_0.3.1       future_1.34.0     foreach_1.5.2     ggdist_3.3.2     
[13] janitor_2.2.0     iarm_0.4.3        hexbin_1.28.5     catR_3.17        
[17] glue_1.8.0        ggrepel_0.9.6     patchwork_1.3.0   reshape_0.8.9    
[21] matrixStats_1.4.1 psychotree_0.16-1 psychotools_0.7-4 partykit_1.2-23  
[25] mvtnorm_1.3-2     libcoin_1.0-10    psych_2.4.6.26    mirt_1.43        
[29] lattice_0.22-5    eRm_1.0-6         lubridate_1.9.4   forcats_1.0.0    
[33] stringr_1.5.1     dplyr_1.1.4       purrr_1.0.2       readr_2.1.5      
[37] tidyr_1.3.1       tibble_3.2.1      ggplot2_3.5.1     tidyverse_2.0.0  
[41] kableExtra_1.4.0  formattable_0.2.1

loaded via a namespace (and not attached):
  [1] splines_4.4.1        bitops_1.0-6         R.oo_1.27.0         
  [4] cellranger_1.1.0     deSolve_1.40         rpart_4.1.23        
  [7] lifecycle_1.0.4      rprojroot_2.0.4      globals_0.16.3      
 [10] MASS_7.3-61          backports_1.5.0      magrittr_2.0.3      
 [13] vcd_1.4-13           Hmisc_5.2-1          rmarkdown_2.29      
 [16] yaml_2.3.10          sm_2.2-6.0           sessioninfo_1.2.2   
 [19] pbapply_1.7-2        RColorBrewer_1.1-3   abind_1.4-8         
 [22] audio_0.1-11         expm_1.0-0           quadprog_1.5-8      
 [25] R.utils_2.12.3       RCurl_1.98-1.16      pracma_2.4.4        
 [28] nnet_7.3-19          listenv_0.9.1        testthat_3.2.2      
 [31] RPushbullet_0.3.4    vegan_2.6-8          parallelly_1.40.1   
 [34] svglite_2.1.3        permute_0.9-7        codetools_0.2-19    
 [37] xml2_1.3.6           tidyselect_1.2.1     farver_2.1.2        
 [40] base64enc_0.1-3      jsonlite_1.8.9       polycor_0.8-1       
 [43] ks_1.14.3            progressr_0.15.1     Formula_1.2-5       
 [46] survival_3.7-0       systemfonts_1.1.0    tools_4.4.1         
 [49] gnm_1.1-5            ragg_1.3.3           snow_0.4-4          
 [52] Rcpp_1.0.13-1        mnormt_2.1.1         admisc_0.36         
 [55] xfun_0.49            here_1.0.1           mgcv_1.9-1          
 [58] msm_1.8.2            distributional_0.5.0 ca_0.71.1           
 [61] withr_3.0.2          beepr_2.0            fastmap_1.2.0       
 [64] fansi_1.0.6          digest_0.6.37        timechange_0.3.0    
 [67] R6_2.5.1             textshaping_0.4.1    colorspace_2.1-1    
 [70] R.methodsS3_1.8.2    inum_1.0-5           utf8_1.2.4          
 [73] generics_0.1.3       renv_1.0.10          data.table_1.16.4   
 [76] SimDesign_2.17.1     htmlwidgets_1.6.4    pkgconfig_2.0.3     
 [79] gtable_0.3.6         lmtest_0.9-40        pcaPP_2.0-5         
 [82] brio_1.1.5           htmltools_0.5.8.1    scales_1.3.0        
 [85] snakecase_0.11.1     irtoys_0.2.2         rstudioapi_0.17.1   
 [88] tzdb_0.4.0           checkmate_2.3.2      nlme_3.1-165        
 [91] curl_6.0.1           zoo_1.8-12           KernSmooth_2.23-24  
 [94] relimp_1.0-5         vcdExtra_0.8-5       fda_6.2.0           
 [97] foreign_0.8-86       pillar_1.9.0         vctrs_0.6.5         
[100] Deriv_4.1.6          cluster_2.1.6        dcurver_0.9.2       
[103] GPArotation_2024.3-1 htmlTable_2.4.3      evaluate_1.0.1      
[106] cli_3.6.3            compiler_4.4.1       rlang_1.1.4         
[109] future.apply_1.11.3  labeling_0.4.3       mclust_6.1.1        
[112] fds_1.8              plyr_1.8.9           rainbow_3.8         
[115] stringi_1.8.4        hdrcde_3.4           viridisLite_0.4.2   
[118] munsell_0.5.1        Matrix_1.6-5         qvcalc_1.0.3        
[121] hms_1.1.3            ltm_1.2-0            PerFit_1.4.6        
[124] readxl_1.4.3        

15 References